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Abstract. In this paper we present a hybrid approach to numerically solve two-dimensional 
electromagnetic inverse scattering problems, whereby the unknown scatterer is hosted by a possibly 
inhomogeneous background. The approach is 'hybrid' in that it merges a qualitative and a quantita- 
tive method to optimize the way of exploiting the a priori information on the background within the 
inversion procedure, thus improving the quality of the reconstruction and reducing the data amount 
necessary for a satisfactory result. In the qualitative step, this a priori knowledge is utilized to im- 
plement the linear sampling method in its near-field formulation for an inhomogeneous background, 
in order to identify the region where the scatterer is located. On the other hand, the same a priori 
information is also encoded in the quantitative step by extending and applying the contrast source 
inversion method to what we call the 'inhomogeneous Lippmann-Schwinger equation': the latter 
is a generalization of the classical Lippmann-Schwinger equation to the case of an inhomogeneous 
background, and in our paper is deduced from the differential formulation of the direct scattering 
problem to provide the reconstruction algorithm with an appropriate theoretical basis. Then, the 
point values of the refractive index are computed only in the region identified by the linear sampling 
method at the previous step. The effectiveness of this hybrid approach is supported by numerical 
simulations presented at the end of the paper. 
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1. Introduction. Acoustic, elastic or electromagnetic scattering is a physical 
phenomenon where an incident wave is scattered by an inhomogeneity (or an obsta- 
cle) and the total field at any point is represented by the sum of the incident and 
the scattered field. From a mathematical viewpoint, in the direct problem the phys- 
ical parameters and the geometry of the inhomogeneity are given and the unknown 
is represented by the scattered field. In the inverse scattering problem, one aims at 
recovering the shape and the physical properties of the object of interest from mea- 
surements of the scattered field. The direct scattering problem is in general well-posed 
in the sense of Hadamard and therefore its approximate solution can be determined by 
means of stable numerical methods. On the contrary, the inverse scattering problem 
is ill-posed in the sense of Hadamard (specifically, the unknown physical parameters 
of the scatterer are mapped onto the measured scattered field by a compact operator), 
and its solution must be addressed by means of some regularizing approach [TT] . 

Most inverse scattering methods belong to two different sets of algorithms: the 
family of qualitative approaches and the family of quantitative approaches. Qualita- 
tive methods [8] require the regularized solution of a linear integral equation of the first 
kind, parameterized over a grid of sampling points covering the investigation domain. 
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The Euclidean norm of such regularized solution behaves as an indicator function 
of the unknown scatterer, since it is bounded when the sampling point is inside the 
target, grows up if the point approaches its edge and can be made arbitrarily large 
when it is taken outside. The advantages of qualitative methods are that they are fast 
and need very few a priori information to work. However, the informative content of 
their output is rather poor, since they are essentially visualization techniques and do 
not provide quantitative reconstructions of the inhomogeneity. 

On the other hand, quantitative inverse scattering methods [3 [11] are, in general, 
iterative schemes that minimize an appropriate functional, starting from an initial- 
ization mask for the refractive index of the scatterer, and then reconstruct its point 
values after an optimized number of iterations. In principle, they can provide all the 
required information on the problem, although by means of a notable computational 
effort. However, they can suffer from local minima problems, since the number of 
unknowns is typically larger than the data amount at disposal in scattering experi- 
ments, or because often the initialization of the method is not precise enough for a 
proper convergence to the solution. The effectiveness of a quantitative method can 
be notably increased by incorporating some a priori information (when available) on 
the scatterer into the minimization process. Although the most traditional approach 
is to use this knowledge for a better initialization of the algorithm, an improvement 
in its performance can be achieved when such information is integrated into each step 
of the iterative scheme, e.g., by decreasing the complexity of the target [14], or by 
adding appropriate constraints to the minimization technique [28] . 

Recent developments in inverse scattering are concerned with the formulation 
of hybrid methods, i.e., methods merging different techniques in order to integrate 
and optimize the different kinds of information they can provide (see e.g. [6] [7] [18]). 
In particular, a priori knowledge on the scatterer and/or qualitative techniques are 
typically utilized to improve the performance of some quantitative method, e.g., by 
reducing the number of unknowns, or by providing a more precise initialization. In this 
paper, we rely on a rather general strategy for the formulation of hybrid techniques, 
based on the following three steps: 

1. the a priori information on the physical characteristics of the (possibly inho- 
mogeneous) background is coded into the corresponding Green's function; 

2. a segmentation between the background and the scatterer is realized by apply- 
ing a qualitative method, thus reducing the size of the investigation domain, 
i.e., the number of unknowns in the following step; 

3. a quantitative method is applied only in the region highlighted by the quali- 
tative method, in order to reconstruct the point values of the refractive index 
in this region. 

Although this scheme is of general applicability, in this paper we present a specific 
implementation: more precisely, the Green's function of the (inhomogeneous) back- 
ground is numerically computed by means of the method of moments [I5j [25] , which 
solves the forward scattering problem. Then the linear sampling method [13] is applied 
in its near-field formulation for inhomogeneous backgrounds [9j [12] in order to visual- 
ize the scatterer of interest, and a procedure based on active contours [TJ [10] is used 
for extracting the (approximate) shape of its support. Finally, the inverse scattering 
problem is solved by generalizing and applying the contrast source inversion method 
P~4j[27][28] to what we call the 'inhomogeneous Lippmann-Schwinger equation', since 
it is a generalization of the integral formulation of the direct scattering problem to 
the case of an inhomogeneous background: such equation is deduced in our paper as 
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a preliminary but essential tool that allows encoding the a priori information on the 
background into each step of the iterative reconstruction procedure. Note that the 
initial qualitative approach allows computing the point values of the refractive index 
only in the region identified by the linear sampling method, thus reducing the number 
of unknowns to be determined by the quantitative method. 

The plan of the paper is as follows. In Section [2] we introduce the differential 
formulation of the direct scattering problem and recall some known results. Sec- 
tion [3] addresses the same problem from an integral perspective, and as a result the 
inhomogeneous Lippmann-Schwinger equation is derived from the direct differential 
formulation. In Section [4] we recall the key ideas of the linear sampling method in the 
near-field case, and adapt the contrast source inversion method to the inhomogeneous 
Lippmann-Schwinger equation. In Section [5] some numerical examples are described, 
also including an application to breast cancer detection by using microwaves. Finally, 
our conclusions are offered in Section [6l 



2. The scattering problem: differential formulation. We consider a rather 
general two-dimensional and time-harmonic scattering problem, whereby the back- 
ground medium is inhomogeneous and described by a piecewise continuously differ- 
entiable refractive index ITT1 [12] 



n b (x) — 
£0 



/ \ .o-(x) 
e(x) +1^^ 



(2.1) 



In (2.1), x = (xi,X2) is a generic point of R 2 , So > is the dielectric permittivity of 
vacuum, e(x) > 1 and <j(x) > are the point values of the dielectric permittivity and 
conductivity of the medium, uj is the angular frequency, and i = More precisely 

(see e.g. fig. |2.1[ ), we assume that there exists a finite number N + 1 of open and 
connected C 2 -domains ^ Cl 2 , with i = 0, . . . , N, such that 1) ^ n flj = for i j; 
2) R 2 = iCo^; 3) Qi is bounded for each 0; 4) n b \ Qi e C 1 ^) for alH = 0, . . . , TV; 
5) nb(x) = no £ C for x G £lo? with Im{no} > 0; 6) there exists a subset Jo of the 
finite set {1,2,..., N} such that flo := int {x G R 2 : rib(x) = no} = ^0 U (Uj G j Qj) 
(in particular, Jo = if and only if ^0 = ^o)- Finally, we assume that the magnetic 
permeability is constant on all R 2 . Note that a) the domains Qi do not need to 
be simply connected: e.g., in fig. |2.1[ ^3 has ^4 as its hole and, in turn, ^1 has 
holes corresponding to ^2 and Q3 U ^4; b) the boundaries dVti may be either curves 
where a discontinuity of n b occurs, or boundaries of virtual domains Qi immersed in 
a homogeneous region, where to host a scatterer that will be introduced later: this 
trick allows some notational simplifications. 
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Figure 2.1. Scheme of the reference background for the scattering problem. In this case, 

Next, we define the Green's function G(x; y) of the background as the radiating 
solution of the equation [12] 

A x G(x;y) + k 2 n b (x)G(x;y) = -S(x - y) for xGR 2 , (2.2) 

where k = uo/c is the wave number in vacuum, c being the speed of light in free space. 
The existence and uniqueness of G(x;y) (at least for y G f^o) can be proved as in 
[TTJ [12]: in fact, if we denote the Green's function of the homogeneous medium in 
fto by $>(x; y) := \ko \x — y\), where Hq 1 ^ is the Hankel function of first kind of 
order zero, fco := &V^o an d lm{^/no} > 0, we can write 

G(x;y) = *(x;y)+ul(x;y), (2.3) 

where ul(x;y) is the perturbation to $>(x;y) due to the inhomogeneous medium in 
R 2 \ £l . 

We point out that [12] G(-;y) G C 1 (IR 2 \ {y}), i.e., the discontinuities of across 
some of the boundaries dQi only affect the smoothness of second (and higher order) 



derivatives on the boundaries themselves. Actually, eq. (2.2) is to be understood 
as a set of N + 1 equations, one for each domain Q^, linked by proper transmission 
conditions at the boundaries dVti. These conditions are imposed by physics, which 
states the continuity of the tangential components of both the total electric field E 
and the magnetic field H: this is equivalent to the continuity of the vector fields v x E 
and vx H, where ' x ' denotes the vector product and v is the unit normal at a point of 
dVt{. In our two-dimensional setting, Cartesian axes are chosen so that G is the non- 
zero component of the total electric field E = (0, 0, G), which vibrates perpendicularly 
to the scattering plane (i.e., the electromagnetic field is assumed to be TM-polarized). 
Moreover, from the time-harmonic Maxwell equation [11] curl£? — \kH = 0, it follows 
that H = T^(d 2 G, — <9iG, 0), having denoted by dj the partial derivative with respect 
to the variable x^ for j = 1,2. Finally, in the same reference system, the normal 
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v can be written as (^1,^2,0): accordingly, v x E = G, -^1 G, 0) and v x H = 
(0, 0, —v\ d\G — diG). Then, the continuity of v x E and v x H corresponds to the 
continuity of G and its normal derivative v • VG = d v G. 

In our scattering problem, the scatterer is assumed to take up the spatial region 
D = yjfLi^i-, with 1 < M < N; accordingly, the whole propagation medium is 
described by a refractive index n(x) such that, in general, n(x) ^ n^x) for xGdj, 
with i = 1, . . . , M. In any case, we still require that n\^ Li G C 1 ^) for alH = 0, . . . , N. 

Moreover, for a unit point source placed at xq G Clo \ (^0 HD), we denote by u(x; xq) 
the non-zero component of the total electric field E = (0,0,^), which, as before, 
vibrates perpendicularly to the scattering plane. Then, the differential form of the 
scattering problem we are interested in is [12] 



A x u(x; xq) + k 2 n(x) u(x; xq) = —S(x — xq) 
u(x; xq) = u l (x; Xq) + u s (x; xq) 

= 0, 



for x G M 2 



lim 

r— )>oo 



r ( du ° ■ r ■ 

V ~dr ~ 1 



(a) 
(b) 

(c) 



(2.4) 



where u % (x;xq) — G(x;xq) is the incident field, u s (x;xq) is the scattered field and 
( |2.4| ) (c) is the Sommerfeld radiation condition, which holds uniformly in all directions 
x = x/\x\. If Im{no} > 0, such a condition can be relaxed [12 by only requiring the 
boundedness of u s (-;xq) in R 2 . 

Note that, from eq. (2.3) and the identification u l (x]Xo) = G(x;xo), we can 
rewrite eq. (2.4)(b) as u(x;xo) = &(x;xo) + [ul(x;xo) + u s (x;xo)]: then, as be fore , 
the existence and uniqueness of a solution u(-;xq) G C 1 (M? \ {^o}) to problem (2.4) 
can be proved as in [12 , and eq. ( |2.4| )(a) is again to be regarded as a set of N + 1 
equations, one for each domain f^, with the above transmission conditions at the 
boundaries <9£V 

Moreover, since both the incident field u 1 {-]Xq) and the tota l field u(-;xo) are in 
C 1 ^ 2 \{xo}), also the scattered field u s (-; xo), expressed by (2.4)(b) as the difference 
u{-\ xq) —u l {-\ xq) : is in C 1 (R 2 \ {xo}). Actually, it is not difficult to establish further 
regularity properties of u s (-;xo), i.e., u s (-;xo) G C 1 (M 2 ) and Au s (-;xo) G C 1 (f^) for 
all i = 0, . . . , N. To this end, we observe that a simple algebraic manipulation of eq. 



(2.4) (a) yields 



A x u(x; xq) + k 2 nt)(x)u(x; xq) = —S(x — xq) + k 2 fiQ- K 



n(x) 



n 



u(x\x ), (2.5) 



i.e., by setting m(x) := [n&(x) — n(x)]/ho and remembering that fco = &V^o? 

A x u(x; xq) + k 2 nb(x)u(x; xq) = —S(x — xq) + k^m{x)u{x] xq). (2.6) 



By using equations (2.2) and (2.4)(b), as well as the above identification u % {x]Xq) 
G(x]Xq), eq. (|2.6[) is easily seen to become 



A x u s (x; xq) + k 2 nb(x)u s (x; xq) = — —klm(x)u(x; xq) 



(2.7) 



The structure of eq. ( |2.7[ ), which parallels that of eq. ( |2.2[ ), allows regarding u s (-; xq) 
as the third component of the electric field E s = (0, 0, u s ) radiated by the equivalent 
source — kQm(x)u(x; xq) in the background medium described by the refractive index 
rib. We now note that this equivalent source is compactly supported, since suppm 
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coincides with the support D of the scatterer, which is compact by assumption. On 
the other hand, by hypothesis, we have that xq G \ (&o HD), i.e., the point xo 
is at a finite distance from the equivalent source: accordingly, there exists an open 
neighbourhood U XQ C of such that 1) U XQ flsuppm = 0; 2) both the electric field 
E s and the corresponding magnetic field H s = t^(9 2 ^ s , — d\u s , 0) are well-defined and 
bounded in U XQ . Even more, from the transmission conditions imposed by physics 
and previously recalled, we can conclude that both u s (-;xq) and d v u s (-]Xo) are in 
C°(U Xo ). But here the unit normal is arbitrary, since no physical interface, i.e., no 
discontinuity of the refractive index occurs inside U XQ , where in fact rib(x) = fiQ\ 
then, in particular, djU s (-;xo) G C°(U XQ ) for j = 1,2. Accordingly, we have that 
u s (-;xq) G C 1 (/7 ;Eo ). On the other hand, we know that u s (-;xq) G C 1 (R 2 \ {xq}): we 
then conclude that u s (-;x ) G C 1 (R 2 ). 
Finally, from eq. (|2.7|) we find 



A x u s (x; xo) — —k 2 ni ) (x)u s (x; xq) + klm(x)u(x; xq). (2.£ 
Since, by hypothesis, n b ,m G C 1 ^) for all i = 0, . . . , N, while u s (-;x ) G C 1 (R 2 ) 



and u(-;x ) G C X (R 2 \ {x }), from (2.8) we have Au s (-]Xq) G C 1 ^): indeed, the 
singularity of u(x;xo) for x = xo is cancelled out by m(x), since xo ^ suppm by 
assumption, while u s (-]Xo) is bounded at infinity (actually, it is bounded in all R 2 

na). 

Of course, as a particular case (i.e., n(x) = n^x) Vx G R 2 ) of the above discussion, 
the same regularity properties also hold for the field (•;?/) introduced in (2.3), i.e., 



l(-;y) G C X (R 2 ) and Aug (■;!/) G C 1 (^) for alH = 0, . . . , N. 

3. The inhomogeneous Lippmann-Schwinger equation. The goal of the 
present section is to derive an integral equation for the scattered field u s (-;xq) from 



of eq. (2.7), the natural candidate for an integral formulation is the 'generalized' 



the differential formulation (2.4|) of the scattering problem. Observing the structure 
of eq. (2.7), the natural cand 
Lippmann-Schwinger equation 

u s (x;x ) = -fcg / G(x;y)m(y)u(y;x )dy. (3.1) 



In fact, an exact equivalence between differential and integral formulation is proved 
in [11] for the three-dimensional acoustic and electromagnetic cases under strong 
regularity assumptions on the refractive index: the background is assumed to be 
homogeneous (i.e., rib(x) = ho G C for all x G R 3 ), and n must be continuously 
differentiate on the whole space (however, see [17 for a variational approach notably 
relaxing the latter requirement). In particular, under such assumptions, this proof 
can provide a rigorous justification of the following procedure: compute A x u s (x; xq) 
by interchanging the Laplacian operator with the integral in ( |3.1[ ), then replace 

A x $(xi y) = -kfo(x; y) - S(x - y) (3.2) 



(which is the analogous of eq. \2.2\ in the homogeneous case G(x;y) = <&(x;y)) into 
eq. ( |3.1| ) to obtain 

A X U S (X] Xq) + koU S (x] Xq) = — — kQim(x)u(x] Xq) . (3-3) 



Eq. (3.3) suggests regarding the term in square brackets at its right-hand side as an 
equivalent, although unknown, source radiating the field u s (x; xq) in the homogeneous 
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background, which is consistent with the physical interpretation of (3.1 ). In any case 



problem (2.4). 



by means of eq.s (3.2) and (3.3), it is easily verified that u s , as given by (3.1), solves 



However, as correctly pointed out in [22], even a mere discontinuity of the physical 
parameters at the interface dD between the scatterer and the (homogeneous) back- 
ground suffices to invalidate, in general, the previous procedure: as a consequence, 
even with a homogeneous background, the integral formulation becomes more com- 
plicated (owing to the occurrence of boundary terms on dD) and its equivalence with 
the differential formulation is more difficult to prove. 

To our knowledge, in the case of an inhomogeneous background, results concern- 
ing the equivalence between the differential and integral formulation of a scattering 
problem are n ot av ailable. In this section, we shall limit ours elves to deriving the inte- 
gral equation (3.1 ) (for x G f^o) from the differential problem (2.4): indeed, a thorough 
discussion of the equivalence between (3.1) and (2.4) would be very technical and be- 
yond the framework of this paper. However, it is likely that an exact equivalence 
actually holds, since, as recalled above, the scattered field u s is in C 1 (M 2 ): then, the 



boundary contributions on each dfli that would appear in (3.1) from Green's second 



theo rem [11 appl ied i n each domain Qi cancel out, as detailed in the following Lemma 



3.1 



and Theorem 



3.2 



for the case x G ^o- 

Finally, it is worth noting that, even in two dimensions, the previous argument 
fails [22] whenever the magnetic permeability is discontinuous across dQi. This is 
analogous to what happens in the two or three-dimensional acoustic case [22], where 
the pressure field is continuous across the discontinuities of the fluid density, but not 
continuously different iable. 

Lemma 3.1. Let (with i = Q,...,N), £Iq, n^x), G(x;y) be as above, let 



Br := {x G 



p2 . 



< R} be such that Br D ufL^i, and define fi/v+i 



Moreover, let w G C 1 (Br) be such that Aw G C°(f^) for all i 
the following generalization of Green's formula holds: 



1, 



= B R \ug 1 n i . 

,7V + 1. Then, 



w(x) 



dB E 



dw dG(y;x) 
- {y) G(y; X )-w(y)^- r 



ds(y) + 



(3.4) 



/ [Aw(y) 

JB R 



k 2 n b (y)w(y)] G(y;x) dy 



Vx g B R nsi . 



Proof Consider an arbitrary point x G Qr := Br D Clo: since fl R is open, there 
exist p > and a ball B{x\ p) := {y G M 2 : \y — x\ < p} such that B(x;p) C Qr. 
Moreover, G(-;x) solves a particular case of the differential problem (2.4), with the 
identifications xo = x, u(-;xo) = G(-;x), u l (-;xo) = <£(•;#) and n = n^: accordingly, 
remembering eq. ( |2.3[ ) and the regularity properties stated in the previous section, we 
have that G(-;x) G C 1 (B R \{x}), AG(-;x) G C 1 (ti R \{x}) and AG(-;x) G C 1 ^;) 
for alH G {1, . . . , N} \ Jo, where the index set Jo has been defined in assumption no. 
6) soon below eq. ( |2.1| ). Then, given the regularity properties assumed for w, we can 
apply the usual Green's second theorem [TTJ[T6] in the domain Qr \ B{x\ p), i.e., 



Q R \B(x;p) 

UdB(x;p) L 



[G{y\ x)Aw(y) - w(y)A y G(y; x)} dy 



(3.5) 



/ 



' ,dw(y) dG{y;x) 
G(y;x)— w(y) 



dv 



dv(y) _ 



ds(y), 
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as well as in any other domain 1^, with i G {1, . . . , N} \ Jo, i.e., 



[G(y; x)Aw(y) - w(y)A y G(y; x)] dy 



(3.6) 



,dw(y) dG(y;x) 
G(y;x)— w(y) 



dv 



dv{y) 



ds(y). 



In eq.s (3.5) and (|3.6|) , it is understood that the unit normal is chosen as outward 



with respect to each domain. We now sum eq. (3.5) with all the equations (3.6) 
obtained for i G {1,...,7V} \ Jo: note that, except for 8Br and dB(x;p), all the 
boundary integrals are taken twice, with opposite orientation of the unit normal. 
Since the integrand functions are continuous on the boundaries, these integrals cancel 



y 7^ x into (3.5) and (3.6). As a result, we find 



out. Moreover, by eq. (2.2), we can substitute A y G(y;x) = —k ri} ) (y)G(y;x) for 



L 



B R \B(x;p) 



[Aw(y) + k 2 n b (y)w(y)] G(y; x)dy = 



(3.7) 



dB R UdB(x;p) 



,dw(y) dG(y;x) 
G(y;x)— w(y) 



dv 



dv(y) 



ds(y). 



We now focus on the integral over dB(x;p), say IdB{x\p)i a t the right-hand side of 
(3.7): remembering eq. (2.3), we have 



1 dB(x;p) 



dB(x;p) 



,dw(y) d$(y;x 
®(y m ,x) — w(y) 



dB(x;p) 



u s h (y;x 



dv 
dw(y) 
dv 



w(y) 



dv(y) 
du s b (y;x) 



ds(y) + 



(3.8) 



dv(y) 



ds(y). 



Now, the second integral in (3.8) vanishes as p —> 0, since the integrand is bounded 



and the measure of the integration domain tends to zero. As regards the first integral 



in (3.8), we recall [8 the following asymptotic behaviors 



1 



In- 



1 



2tt \x-y\ 
d$(y;x) 1 1 



dv(y) 



2tt \x - y\ 



+ 0(1) 

0(\x-y\\n\x-y\) 



as \x-y\ -^0, (3.9) 
as \x-y\ -^0. (3.10) 



By using ( |3.9| ) and ( |3.10| ) , the integral mean value theorem applied to the first integral 

-w(x) as p — >• 0. Then, thesis (3.4) is 
(|3.7[): indeed, remembering eq. 



in (3.8) easily shows that the latter tends to 
obtained by taking p 



in eq. ( |3.7[ ): indeed, remembering eq. ( |2.3[ ) and the 
regularity of ul(y; x), the singularity of G(y; x) for y — » x is only due to $>(y; x), i.e., 
it is weak and then the integral on Br converges. □ 

Theorem 3.2. Let Qi (with i = 0, . . . ,N + l), VLq, D, n h (x), n{x), m{x), G(x; y) 

be as above, and let xo G flo \ (flo D D) be as above. If u s (-;xo) G C 1 (IR 2 ) ; with 
Au s (-',xo) G C 1 (Clij for all i = 0, ...,iV ; is a solution of the differential problem 
(2.4), then u s (-;xq) solves the integral equation (3.1) for x G ^o- 

Proof. Consider an arbitrary point x G VLq. Let Br := {x G M 2 : \x\ < R} 
be an open disk with exterior unit normal v and radius R large enough, so that 
Br Z> [ufL-^Cli U {x}] . By hypothesis, u s (-]Xq) is regular enough in the domain Br 
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to be represented by means of the generalized Green's formula (3.4) in Qr := BrC\VIo, 
i.e., 



U S (x]Xq) = I 

JdB L 



du s (y;x ) 
dv(y) 



G(y;x) - u s (y;x ) 



dG(y 



dv(y) 



y) J 



ds(y)+ 



(3.11) 



{A y u s (y; x ) + k 2 n b (y)u s (y; x )} G(y; x)dy, x e Q t 



First we prove that the integral on 8Br in (3.11) is zero. To this end, consider 



{x e 



p2 . 



< r} such that r > R, and apply Green's second theorem [TT] in 



B r 

the domain B r \ Br: by choosing the unit normal v as outgoing from both Br and 
B ri and observing that both the fields u s (-;xo) and G(-;x) verify the same Helmholtz 
equation in B r \ Br C Qq (where rib(x) = n{x) = no), we find 



/ 

-I { 

JdB r I 



du s (y;x ) 



dv(y) 
du s (y;x ) 
dv{y) 



G(y;x) - u s (y;x ) 



dG{ y] x) ' 
du(y) . 



ds(y) 



(3.12) 



G(y;x)-u s (y ] x ) d ^^}ds(y). 



We now recall [TT] that any radiating solution v of the Helmholtz equation (with 
generic wave number k > 0) outside a disk in R 2 has the following asymptotic behavior 



v(x) 



Voo(x) + O 



oo, 



(3.13) 



where Voo is the far-field pattern of v. If v 1 and v 2 are two such solutions, from (3.13) 
we find [TT] 



v 1 (x) 



dv 2 (x) 
dr 



n 2\kr 



ik- 



r — » oo, 



(3.14) 



uniformly for all directions. By applying (3.14) to the radiating fields u s (-;xq) and 
G(-;x), we easily find that the integrand function at the right-hand side of ( |3.12 ) is 
O (^2 ) and then the integral itself vanishes as r — >• oo. This is even more true if the 
wave number k > is replaced by k\fn$ with Im { v^o} > 0, since the attenuation of 
the fields and their derivatives is faster. 

As regards the integral on Br in ( |3.11[ ), we come back to eq. (2.7) in Section [2] 
and recall that the singularity of u(x;xo) for x = xo is cancelled out by m(x), since 
xo £ suppm: then, substituting ( |2.7[ ) into ( |3.11 ) and taking the limit as R — >• oo, we 
can write 



u s (x;x ) 



x)m(y)u(y\ x ) dy \/x e Cl , 



(3.15) 



which is immediately written in the form (3.1) by using the reciprocity property 
0H G(y;x) = G(x]y). □ 



We call eq. (3.1) the 'inhomogeneous Lippmann-Schwinger equation', to empha- 



size that the reference background is (or may be) inhomogeneous. In the next section, 
we shall apply an inversion algorithm to this equation in order to compute the point 
values of m, i.e., of the refractive index inside the region under investigation. 
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4. A hybrid scheme. In general, iterative methods for the quantitative solution 
of inverse scattering problems are applied to the homogeneous Lippmann-Schwinger 
equation. Such methods take as input the scattering data collected by antennas placed 
outside the investigated area and inside a homogeneous zone, and provide as output 
the reconstruction of the refractive index everywhere in the inhomogeneous region. 
The main drawback of this computational approach is that, particularly when the 
scattering experiment is performed with a single and fixed frequency, the number of 
unknowns is typically much larger than the number of measured data, and therefore 
the reconstruction accuracy is often rather low. However, there are situations where 
just a certain part of the inhomogeneous region is unknown and of interest for practical 
applications, while information is available about the refractive index of the rest of the 
domain. In this case, the quantitative inverse scattering method can be applied to the 
inhomogeneous Lippmann-Schwinger equation, provided one is able to compute the 
Green's function of the inhomogeneous (and known) background and to approximately 
identify the region taken up by the scatterer under investigation. In the following of 
the current section, we describe an implementation of this approach essentially based 
on the contrast source inversion method. Therefore, in order to realize the proposed 
hybrid approach, the ingredients we need are: 1) a method for computing the Green's 
function of the background; 2) a qualitative method to visualize (an overestimate of) 
the scatterer support inside the background itself; 3) a quantitative scheme for the 
inversion of the inhomogeneous Lippmann-Schwinger equation in the region identified 
at step 2). 

4.1. The computation of the Green's function. A handy way to compute 
the Green's function of an inhomogeneous background is given by the method of 
moments (MOM) [15] [25]. Since G(-;xo) solves a particular case of problem ( |2.4[ ) 
with the identifications u(-;xq) = G(-;xo), u z (-;xo) = <£(•; xo) and n(x) — rib(x), it 
also satisfies the homogeneous Lippmann-Schwinger equation [TTJ [12] 



G(x;x ) = <5>(x;x ) -kl / <S>(x;y)m b (y)G(y; x ) dy for x G M 2 , (4.1) 

with rrib(x) := [ho — rib (x)]/ ho and xo G Since suppm^ is compact, the integration 
domain in ( |4.1[ ) is bounded: then, we can consider a finite (and not necessarily uni- 
form) partition of supp by L cells Ai, . . . , Al, chosen so small that and G(-]Xq) 
can be assumed to be constant inside each cell. Then, eq. ( |4.1[ ) can be approximated 

as 

G(x;x ) ~ $(x;xo) - kl V m 6 (%)G(^; x ) / <$>(x;y)dy, (4.2) 

3=1 Ja > 

where yj G M 2 identifies the center of the cell Aj. Furthermore, following [25 , we 
approximate each square cell as a circular cell of the same area, so that the integral 
in eq. (4.2) can be evaluated as 

/ $(x;y)dy= 1 -[ {h\x - y\)dy « (4.3) 

Ja 3 4 J A . 

( rnkpaj ~ (1) 



Ji{koa j )H^ ) {ko\x-y j \) for x £ Aj 



i 

2 L 



nkocijH^ \kodj) + 2i 



for x e Aj, 
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where aj := yj Aj/tt, Aj is the area of the cell Aj, J\ is the Bessel function of first 

order and is the Hankel function of first kind of order one. We point out that the 
method can be notably speeded up by utilizing the fast Fourier transform algorithm 
as it is described for example in [4j [29]: we shall adopt this version of the MOM 
method for our numerical simulations. 

4.2. The linear sampling method. The linear sampling method (LSM) is 
the earliest and most used qualitative method in inverse scattering: in the case of 
an object embedded in a homogeneous background, it provides a reconstruction of 
its support by only knowing the field measured around it. When the background is 
inhomogeneous and its Green's function is known, the LSM can be extended to the 
case of an inhomogeneity immersed in a medium with piecewise constant refractive 
index. The basic idea is to write and to approximately solve the modified far-field 
equation [12] 



{u s (x; x ) - u s b (x] x )]g z (x ) dx = G(x; z) for x G I\ (4.4) 

where T := {x e M? : \x\ = Rr} C Ho is the circle of radius Rr where emitting and 
receiving antennas are placed, and z G M 2 is a sampling point inside the investigation 



domain enclosed by T. In (4.4), the Green's function G(x;z) and the field ul(x;xo) 
can be computed by exploiting the knowledge of the background and applying the 
MOM method, while u s (x;xo) represents the measurements, i.e., the data of the 
inverse scattering problem. 

In fact, it can be proved [9j [12] that there exists an approximate solution of 



(4.4) whose L 2 (T)-norm is bounded inside the scatterer, grows up as z approaches 
its boundary from inside and remains very large outside, thus behaving as an indi- 
cator function for the support D of the scatterer itself. This result inspires a simple 
algorithm that approximately solves the modified far-field equation by means of a 
regularization procedure. More precisely, the LSM requires the choice of a numerical 
grid covering the region where the scatterer is placed; then, for each point z of the 
grid, a discretized version of the modified far-field equation is solved by using the 
Tikhonov regularization method [26], and the Euclidean norm of the regularized so- 
lution is plotted for each z. As a result, the boundary of the scatterer is highlighted 
by the points of the grid corresponding to the largest increase of the Euclidean norm. 
We recall that the computational time of this algorithm can be notably reduced by 
applying a no-sampling formulation [2], called 'no-sampling linear sampling method' 
(NSLSM), which is the one adopted in the present paper. 

Finally, we remark that, for the purposes of our implementation, a univocal iden- 
tification of the shape of the scatterer is performed by post-processing the NSLSM 
through an active contour technique (TJ [10] : for reasons that will be clear in the next 
subsection, we choose the parameters of this edge detection algorithm in such a way 
that a slight overestimate of the scatterer is favoured. 

4.3. The contrast source inversion method. In our hybrid scheme, the con- 
trast source inversion (CSI) method [14] [27, 28 is the technique we use for the quanti- 
tative inversion of the inhomogeneous Lippmann-Schwinger equation. More precisely, 
the CSI method computes the point values of the refractive index n(x) inside an in- 
vestigation domain T containing the scatterer support D = supp m. The idea at the 
basis of the CSI method is to split an inverse scattering problem described by an 



integral equation such as (3.1) into two different problems, one defined inside T and 
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the other one on the curve T where antennas are placed. More formally, remembering 
( |2.4| )(b), equation (3.1) is equivalently recast in the form 



and then replaced by the system 

j w(x] xq) = m(x)u l (x] xq) — m(x) [Q T w(-; xo)] (x) for x G T (a) 

\ u s (x;x ) = — [G F w(-;xo)] (x) for x G T, (b) 

where k;(x;xq) := m{x)u(x\ xq), the operator Q T : L 2 (T) L 2 (T) is defined as 



(4.6) 



[£ T /] (x) := %l G(x; z)f(z) dz Vx € T, 



and G r : L 2 (T) -»• L 2 (r) is defined as 



^ / G(x; z)/(,z) dz \/x e T. 



(4.7) 



(4.; 



Note that eq. ( |4.6[)(a) is obtained from eq. (4.5) by multiplying both members for 
m(x), while eq! (|4.6|)(b) is the restriction of eq. (3.1) for x G Y: in particular, 



u s (x] xo) for x G T represents the data of the inverse scattering problem. 
The CSI method consists of minimizing the functional 



F(w, m) 



\\u s -g v w\\ 2 T \\mij 



m G T w\\?n 



\\w\\l 



(4.9) 



where we dropped the dependence on x, xq, and || • ||r, || • \\t denote the L 2 -norms 
on the spaces L 2 (T), L 2 (T) respectively. The minimization of F is performed by 
applying gradient algorithms [24 and the updates are alternately computed for m 
and w [271 EH] • The information on the background is again coded into the Green's 
function, i.e., into the operators Q T and Q v . 

In the standard case of a homogeneous medium, the implementation of the ap- 
proach described above is well-established and does not require further discussion. 
Instead, in the inhomogeneous case two issues need to be addressed. First, we ob- 



serve that eq. (4.6) (a) is written for x G T; on the other hand, the integral equation 
(3.1), whence eq. (|4.6|)(a) should derive, has been proved in Theorem 3.2 only for 



x G an d, in general, T is not contained in Qq. To overcome this drawback, we 
observe that the purpose of the CSI algorithm is to compute the values of n(x) for 
x G T, and these values do not depend on the reference background inside T. Then, 
we can think that the original background in T is replaced by an artificial one, which 
hosts the same homogeneous medium occupying Qq. In other words, this amounts 
to replacing the refractive index n& with hb : such that hb(x) = ho for all x G T and 
hh(x) — r ih(x) for a ll x £ T. Of course, the background Green's function G{x\ z) in 



eqs. (4.7) and (4.8) must be replaced accordingly, i.e., by G{x\z). We point out that 
a good choice of the investigation domain T is the region highlighted by the NSLSM, 
as explained in Subsection |4.2| indeed, such region is an overestimate of D, i.e., it 



contains D, but is also close to be as small as possible, thus minimizing the number of 
points where n(x) is to be computed, i.e., the number of unknowns of the problem. As 
a result, also the number of measurements necessary for a successful implementation 
of the method is optimized. 
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This same trick, i.e., changing the background as just explained, is also useful 
to address the second issue, which is concerned with the computation of the Green's 



function G(x;z). Indeed, eq. (4.7) shows that both x and z vary in T: then, in 



particular, also the singular case x = z is of interest. Now, G{x; z) satisfies eq. (4.1) 



(with the identification z = xq), which implies that if z G suppm^ and x = z, the 



integrand function at the right-hand side of (4.1) is affected by a double singularity 



for y = z: one due to the other one due to G(y;z). This prevents us from 



applying the usual approximation scheme outlined in eq.s ( |4.2[ ) and (4.3): accordingly, 
an ad hoc quadrature rule should be created to numerically compute the integral. 
Such problem is avoided by the previous choice of the new background, i.e., of 



indeed, the region T is now erased by the actual integration domain in (4.1), since 
fhb(x) := [ho — hb(x)]/ho vanishes inside T, while z G T. 

We are now ready to present in the next section some numerical results obtained 
by implementing the theoretical framework developed so far. 

5. Numerical examples. The aim of this section is twofold. First, we numeri- 
cally validate the CSI algorithm in the inhomogeneous case: more precisely, we com- 
pare the performance of the CSI based on the inhomogeneous Lippmann-Schwinger 
equation with that of the traditional CSI (the two algorithms will be called 'inhomo- 
geneous CSI' and 'homogeneous CSI', respectively). To this end, we perform a pre- 
liminary set of simulations by using simple numerical phantoms and without adding 
noise to measurements. Second, we test the effectiveness of the whole hybrid scheme 
presented above, as well as its worthwhileness in real applications, by implementing 
it with a realistic phantom of a female breast slice: the goal is to highlight the pres- 
ence of a tumoral mass inside the inhomogeneous background formed by the healthy 
tissues. 

5.1. Homogeneous and inhomogeneous CSI: a comparison. The first 
set of numerical simulations is concerned with the reconstruction of the two (non- 



absorbing) phantoms visualized in panels (a) and (b) of Figure 5.1 In particular, we 
are going to show an example where the homogeneous CSI fails to provide a satisfac- 
tory result, while the inhomogeneous CSI succeeds: this is obtained by coding some 
additional information into the background, which then becomes inhomogeneous. 



Indeed, consider first the phantom in Figure 5.1 'a): the pixel values of the rela- 
tive dielectric permittivity e r and the geometry of the objects are plotted in Figure 
5.2[ a). The scattering experiment is simulated by choosing a wavelength A equal to 
the length unit adopted for all panels in Figure [5^1 and by using 30 unit point sources 
uniformly placed on a circle centred at the origin and with radius 3 A. The scattered 
field is computed by means of a MOM code, as outlined in Subsection |4.1[ at 30 
points obtained from the previous ones after a rotation of 7r/30. Then, we apply the 
homogeneous CSI (implemented in the error-reducing version of [28] and initialized 
by backpropagation) : as a result, we obtain the satisfactory reconstruction shown in 



Figure 5.2 'b). 



However, if we consider the more complex scenario of Figure 5.1 'b), with pixel 



values of e r plotted in Figure 5.2 'c), the same homogeneous CSI provides the un- 
satisfactory reconstruction shown in Figure (!T2| d): it is evident that the algorithm 
converges to the wrong local minimum. A feasible way to overcome this drawback is 
to consider the square barrier surrounding the two disks as a part of the background: 
this situation is represented in Figure (IT2|e). Then, we implement the inhomogeneous 



CSI: in particular, the reconstruction is now performed only inside the barrier, i.e., 



inside the green square highlighted in Figure 5.2 'f). As a result, the effectiveness of 
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(a) (b) 

Figure 5.1. Visualization of the relative dielectric permittivity e r of the numerical phantoms 
utilized for the first set of scattering experiments. The phantoms are piecewise homogeneous and 
purely dielectric, i.e., non- absorbing. 



the algorithm in reconstructing the two disks is restored, as shown by Figure 5.2 'f) 
itself. 

5.2. The hybrid approach validated against a realistic numerical breast 
phantom. In the following numerical examples we test the effectiveness of the hybrid 
scheme described in Section [4] when applied to scattering situations that might occur 
in practice. To this end, we consider a realistic numerical phantom of a female breast 
slice, whose pixel values of relative permittivity e r and conductivity a are obtained 
from the segmentation of MRI (Magnetic Resonance Imaging) images [20] (such phan- 
toms are downloadable from the web site uwcem.ece.wisc.edu/home.htm). Then, a 
circular tumor with a diameter of 0.5 cm is artificially inserted into the phantom: we 
refer e.g. to [21] for the values of the electric parameters of tumoral tissues at various 
frequencies. Our numerical experiments are performed at a frequency of / = 3 GHz: 



accordingly, in panels (a) and (b) of Figure 5.3 we plot the pixel values of the relative 



permittivity e r and conductivity a of the breast tissues at such frequency. As shown 
in this figure, we also assume that the breast slice is immersed in an infinite and 
homogeneous coupling medium [23]: this trick reduces reflection phenomena, thus 
favouring the penetration of the wave into the breast itself. Finally, we choose 30 unit 
point sources uniformly placed on a circle surrounding the breast, concentric with it 
and having a radius of 6 cm. By using a MOM code, the scattered field is computed 
at 30 points obtained from the previous ones after an angular shift of 7r/30, then these 
field values are perturbed by 3% Gaussian noise. 

A whole and precise reconstruction of such a heterogeneous scatterer as a breast 
slice is extremely difficult to be achieved: the main issue is the strong variability of 
the refractive index from point to point, as well as its large range of values. Then, in 
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order to properly describe the heterogeneity of the tissues, the investigation domain 
should be discretized by using a large number of pixels, which corresponds to a large 
number of unknowns in the CSI algorithm. In practical applications, this number can 
be too large with respect to the amount of data, thus impairing the reliability of the 
reconstruction. Therefore, an inhomogeneous-background approach to the problem 
can be helpful: indeed, when a priori known, a part of the scatterer can be included 
in the background, so that the reconstruction can be performed on a fewer number of 
pixels (i.e., in a smaller investigation domain), which reduces the number of unknowns 
and increases the quality of the reconstruction itself. 

In the first numerical example, whose results are collected and visualized in Figure 
|5.4[ we assume that the known background is formed by the fat internal tissue, the 
skin layer and the coupling medium. In general, we may think that such preliminary 
information on the tissues is obtained from the segmentation of the most recent MRI 
image of the patient's breast, while the properties of the coupling medium are known a 
priori. The pixel values of the relative permittivity e r and conductivity a of these me- 
dia are shown in panels (a) and (b) of Figure [5~4| Since our goal here is to reconstruct 
both the glandular and tumoral tissues inside the breast fat, such background panels 
have been respectively obtained from panels (a) and (b) of Figure |5.3| by removing 
the artificial circular tumor (i.e., by restoring the original healthy tissue) and by ar- 
tificially replacing the central glandular mass with an 'average fat' (i.e., by replacing 
the pixel values of e r and a in the glandular tissue with constant values computed by 
averaging over the pixel values of e r and a in the fat tissue). Then, this background 
is utilized to implement the NSLSM, the qualitative technique outlined in Subsection 
In Figure [53c) we show the visualization of the unknown region as provided by 



4.2 



the NSLSM: more precisely, we plot the values of the indicator function, chosen as 
the reciprocal of the squared norm of the regularized solution to the modified far-field 



equation (4.4). Next, we apply an active contour technique [TJ [10] to this visualiza- 
tion map: as a result, we extract the support of the unknown region T, which is the 
homogeneous inner domain shown in panel (d) of Figure |5.4| More precisely, in this 
panel we plot the pixel values of e r characterizing the artificial background we need to 
consider for applying the inhomogeneous CSI in the investigation domain T (the plot 
of a would be analogous and is not shown here). Indeed, we remind that, in order 
to avoid the theoretical and computational problems described in Subsection |4.3[ the 
domain T should be considered as filled with the outmost medium (i.e., in our case, 
with the coupling medium). The reconstructions provided by the inhomogeneous CSI 
inside T for the relative permittivity and for the conductivity are presented in panels 
(e) and (f) of Figure [5~4| respectively: the results are in rather good agreement with 



the true phantoms of Figure 5.3 



In the second numerical example, we utilize the inhomogeneous CSI to investigate 
the nature of a small inhomogeneity detected by the NSLSM. In this case, while 
maintaining the coupling medium as before, the whole healthy breast slice, just as 
provided by the segmentation of the MRI image, is assumed as background, and the 
reconstruction is performed only on the pixels where the presence of the (prospective) 
tumor is highlighted (for a detailed analysis of the performance of the NSLSM in this 



framework, see [5 ). As before, panels (a), (b) and (c) of Figure 5.5 are concerned with 
the NSLSM: (a) and (b) are the plots of the relative permittivity and conductivity of 
the background, while (c) shows the NSLSM output, whence the investigation domain 
T (containing the inhomogeneity to be analyzed) is extracted by using active contours. 
In panel (d) of Figure [5^5] we plot the relative permittivity of the artificial background 
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(a) (b) 

Figure 5.3. Numerical phantom of the breast slice utilized in our scattering simulations, per- 
formed at a frequency of f = 3 GHz. Panel (a): pixel values of the relative dielectric permittivity 
£ r . Panel (b): pixel values of the conductivity a. In both panels, the tumoral mass is the small disk 
centered at (0.02, 0) m and with a diameter of 0.5 cm. 



utilized for the inhomogeneous CSI: the coupling medium is now inserted into the area 
T detected by the NSLSM. Finally, panels (e) and (f) show the reconstructions of the 
relative permittivity and conductivity provided by the inhomogeneous CSI inside T: 



again, a comparison with the true phantoms of Figure [573] highlights that the quality 
of these reconstructions is rather good. 

6. Conclusions and further developments. In this paper we have proposed a 
hybrid method to numerically solve a large class of two-dimensional inverse scattering 
problems. According to this hybrid approach, a qualitative method, i.e., the NSLSM, 
is first applied to approximately identify the region of primary interest inside an inho- 
mogeneous background; then, a quantitative algorithm, i.e., the inhomogeneous CSI, 
is implemented only in this region to compute the point values of the electric param- 
eters (permittivity and conductivity) of the unknown scatterer. Both steps are made 
possible by the knowledge of the inhomogeneous background: inserting such informa- 
tion in the whole procedure improves the quality of the reconstruction and optimizes 
the number of data necessary for a satisfactory result, as shown by the numerical 
simulations in Section [5] The first (qualitative) step consists of a no-sampling imple- 
mentation of the LSM for an inhomogeneous background, which can be performed 



according to the guidelines recalled in Subsection 4.2 and detailed in the papers cited 
there. Instead, the second (quantitative) step requires a preliminary theoretical in- 
vestigation, in order to extend the usual Lippmann-Schwinger equation to the case of 
an inhomogeneous background. Such a generalization has been obtained in Section [3] 
as a consequence of the differential formulation of the scattering problem: under ap- 
propriate assumptions on the refractive index and on the magnetic permeability, the 
new integral equation maintains its simplest form, i.e., involves no boundary terms 
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across the discontinuities of the electric parameters. Accordingly, the homogeneous 
and inhomogeneous Lippmann-Schwinger equations have the same structure: this has 
suggested the idea of generalizing the CSI algorithm, originally based on the former 
equation, to the case of an inhomogeneous background, handled by the latter. 

As regards possible further developments, a first task is trying to prove the exact 
equivalence between differential and integral formulation of the scattering problem 
(possibly by adopting a variational approach, i.e., by interpreting system (2.4) in a 
weak sense and looking for its solution in an appropriate Sobolev space, cf. [17]): 
as pointed out in Section [3J such equivalence is likely to hold. As a generalization, 
the case of non-constant magnetic permeability ji(x) would also deserve an analogous 
investigation: the corresponding integral formulation would be more complex, owing 
to the presence of boundary terms across the discontinuities of n(x). Of course, the 
same problems can also be formulated in a three-dimensional framework. An ideal 
goal might be to parallel and generalize to the case of an inhomogeneous background 
the results obtained in [22] for a homogeneous one. This would allow extending our 
hybrid approach to very general three-dimensional inverse scattering problems, where 
the issues of a good balance between the amount of data and the number of unknowns, 
as well as the related problem of computational costs, are even more important than 
in two-dimensional set-ups. 
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